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We use advanced statistical tools of time-series analysis to characterize the dynamical complexity 
of the transition to optical wave turbulence in a fibre laser. Ordinal analysis and the horizontal 
visibility graph applied to the experimentally measured laser output intensity reveal the presence of 
temporal correlations during the transition from the laminar to the turbulent lasing regimes. Both 
methods unveil coherent structures with well defined time-scales and strong correlations both, in 
the timing of the laser pulses and in their peak intensities. Our approach is generic and may be 
used in other complex systems that undergo similar transitions involving the generation of extreme 
fluctuations. 
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Fibre lasers are important practical devices that rep¬ 
resent complex nonlinear systems with many degrees of 
freedom [ll-l3|. Typically, the output of a fibre laser in¬ 
volves nonlinear interactions of millions of longitndinal 
cavity modes in regimes far from thermal equilibrium Q . 
In general, wave dynamics in fibre lasers is highly com¬ 
plex, as in other optical wave turbulence systems 
Though the underlying physical effect, nonlinear four- 
wave mixing, is purely deterministic and well understood, 
it is also well known that a deterministic dynamical de¬ 
scription is not fully adequate in this problem, and sta¬ 
tistical tools such as entropy and complexity measures 
should be used to characterize the complex fluctuations 
in the generated output signals. Within this framework 
of wave turbulence, the role of “temperature” is played 
by optical noise that occurs in the gain medium, which 
in fibre lasers leads to an effective “nonlinear noise” due 
to four-wave-mixing. 

Recently, the analogy between hydrodynamic transi¬ 
tion to turbulence and change of operational regimes in 
fibre lasers has been studied both experimentally and 
theoretically Q- A transition from highly ordered las¬ 
ing regime to more irregular lasing, characterized by ex¬ 
treme, apparently random intensity fluctuations was re¬ 
ported. Such transition, being a relevant example of a 
phase transition in a one dimensional physical system, 
was shown to be accompanied by the occurrence of co¬ 
herent spatio-temporal structures Q- 

In this work we address an important question rele¬ 
vant to practical identification of such structures: are 
there underlying correlations and/or specific time-scales 
in the easily measurable intensity fluctuations of laser 
radiation? In order to investigate this issue we use two 
nonlinear analysis tools: ordinal analysis and the hor¬ 
izontal visibility graph 0- We show that both methods 


allow to clearly identify the presence of long-range tem¬ 
poral correlations in the experimentally measured laser 
output intensity. 


In our experiments, we measure an output temporal in¬ 
tensity dynamics of a quasi-CW Raman fiber laser formed 
of 1 km of normal dispersion fiber placed between two 
fiber Bragg gratings acting as cavity mirrors Q- State- 
of-the-art experimental capabilities allowed us to regis¬ 
ter extremely long time traces with total number of in¬ 
tensity data points of 50 million. Taking into account 
the discretisation time of 12.5 ps, the intensity dynamics 
over 625 /rs could be captured. In order to be able to 
compare among time-series recorded at different pump 
power, each time-series is normalized to have zero-mean 
and unit variance. Depending on power, the generation 
regime can be considered as ’’laminar” or ’’turbulent” 
7|. The transition occurs at pump power 0.9 W (see 
7| for details). Despite the radically different coherence 
properties of radiation in these two regimes, the output 
intensity, I{t), looks similar and irregular at all powers, 
as seen in Fig. 1(a), with typical intensity probability 
distribution function (pdf) of intensity values, p{I(t)), 
shown on Fig. 1(b). 


As a starting step of our analysis, we generate two 
new data sets from experimentally measured intensity 
dynamics, I{t). In particular, we filter out from the ini¬ 
tial intensity dynamics only local intensity peaks and cre¬ 
ate a sequence of intensity peak heights, {/max,i}i which 
are above a certain threshold, also shown in Fig. 1(a). 
Because each time series is normalized to zero mean and 
unit variance, the thresholds used are in units of the stan¬ 
dard deviation, cr: with threshold 0, the peaks considered 
are only those above the mean value; with threshold 1, 
the peaks are only those above cr. The number of the 
peaks found in intensity dynamics measured at different 
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FIG. 1. (a) Stochastic dynamics of a quasi-CW Raman fibre 
laser: intensity time-series measured at pump power 0.9 W 
(the dots indicate the peak values above a threshold, indicated 
with a dashed line, and the inset displays a detail), (b) Pdf 
of the intensity values, (c) Schematic representation of the 
two analysis methods: the values {xi} in a time series are 
represented with vertical bars, the ordinal patterns formed 
by {xi,Xi+i,Xi+ 2 ) are indicated in red, the links in horizontal 
visibility graph are indicated with arrows, and the numbers 
indicate the degree (the number of links) of each data point 
(graph node), (d) Number of intensity peaks vs the pump 
power for various thresholds (threshold values are shown by 
different colours). 


pump power level is shown on Fig. 1(d) depending on the 
threshold value. We note that, if the threshold is low, the 
number of peaks decreases with the pump power, but if 
the threshold is high, the number of peaks increases with 
the pump power until 0.9 W, where it is maximum, then, 
with further increase of the pump power, the number of 
peaks diminishes. The transition between two different 
generation regimes — laminar and turbulent regime — 
could be easily detected at 0.9 W. In the following we 
fix the threshold equal to 2, i.e., we analyze only the in¬ 
tensity peaks that are above 2o-. As each local intensity 
peak has also a time instant Ti at which it occurs, we 
generate in parallel another data set: a sequence of time 
intervals between local intensity peaks, {AT^}. 

To investigate the complexity of the dynamics and 
to uncover hidden temporal correlations we use two 
methods of time-series analysis, which are represented 
schematically in Fig. 1(c). The first one, known as or¬ 
dinal analysis Q, transforms a time series {xi} (where 
{xi} is either the sequence of peak heights, {/max.i}, or 
the sequence of time intervals between peaks, {ATi}) into 
a sequence of symbols (referred to as ordinal patterns, 
OPs), by considering the order relation among D values 
of the time-series. For example, there are 2 different or¬ 


dinary patterns of length 2: pattern ‘OF if Xi < Xi+i and 
pattern ‘10’ if xt > Xi+i. If H = 3, there are 6 possible 
patterns: Xi < Xi+i < Xi +2 gives ‘012’, Xi +2 < Xi+i < Xi 
gives ‘210’, etc. The number of possible patterns of the 
given length D is Dl. In this way, a sequence of patterns 
could be generated from the sequence of the peak heights 
or from the sequence of time intervals between the peaks. 

After defining the sequence of patterns, one can calcu¬ 
late the probability to find the given pattern in the data 
set, Pi- The entropy computed from their probabilities. 
Pi, of occurrence in the time series, SpE = —J2Pi logPi, 
known as permutation entropy, has been shown to be 
an appropriated measure of the complexity of a time- 
series ]^, 11 , 1 ^ . When there are no serial correlations in 
the time-series {xi}, then all the patterns are equally 
probable and SpE ~ logH!. On the contrary, when 
there are serial correlations, then the OPs are not all 
equally probable, and the permutation entropy will be 
SpE < logi?!. In the following we refer to the normal¬ 
ized entropy, Spe/ logH!, as PE entropy. Thus, with an 
appropriate choice of pattern length D, the OP proba¬ 
bilities and the PE entropy will capture the existence of 
underlying correlations in the time series. 

To verify independently presence of correlations, we 
also apply the so-called horizontal visibility graph (HVG) 
method [l0|. In this approach a time series, {xi}, is con¬ 
verted into a graph by considering each data point, Xi, 
as a node. Any two nodes are connected by a link (or 
edge) if horizontal visibility exists between them [see Fig. 
1(c)]: Xi and Xj are connected if it is possible to trace a 
horizontal line linking Xi and Xj not intersecting interme¬ 
diate data; mathematically, this means that Xi and Xj are 
connected if: Xi, xj > x„ for all i < n < j. Note that 
this graph representation of the time series {xi\ takes 
into account both, the order and the values of the data 
points. Time series with different dynamics are mapped 
into graphs that exhibit distinct topological structures 
0. The topology of a graph is characterized by the de¬ 
gree distribution, p(k), that is the probability that a node 
has k links. Thus, the entropy of the degree distribution, 
Shvg = —^Pk^ogpk (in the following, referred to as 
HVG entropy), is another measure of the complexity of 
the time series {xi} 0- The HVG method also allows 
to analyze different time-scales by constructing the graph 
not from all the “raw” data points, but from lagged data: 
{Xi, , Xip2r: ■ • -I* 

While both analysis methods share the common fea¬ 
ture of transforming the time series {xi} into a sequence 
of integer numbers, {ki} (in the OP case, ki G [1,T^!] is 
the pattern label; if ‘01’, ki = 1; if ‘10’, ki = 2, etc.; in 
the HVG case, ki G [l,iV — 1] is the degree of Xi, with 
N being the number of data points in the time series), 
they have important differences: while the OP method 
requires the pre-definition of the length of the pattern 
D, and does not take into account the values of the data 
points, the horizontal visibility graph method does not 
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require to pre-define an analysis length, and considers 
both, the order relation and the actual values of the data 
points. 






FIG. 2. Probabilities ordinal patterns of length D = 3 vs 
the pump power calculated from the sequence of (a) intensity 
peaks and (b) time intervals between consecutive peaks. Per¬ 
mutation entropy (c) and HVG-entropy (d) calculated from 
the sequence of intensity peaks. 

Figures 2(a) and 2(b) display the probabilities of the 
six patterns of length 3 calculated from the sequence of 
intensity peaks and from the sequence of the time inter¬ 
vals, respectively. We observe that the variation of the 
probabilities with the pump power captures the transi¬ 
tion between two dynamical regimes: below the tran¬ 
sition the OPs are equally probable, while during the 
transition from laminar to turbulent regime their prob¬ 
abilities are different from equiprobability. It can also 
be noticed that, during the transition, the patterns 012 
and 210 become more probable (less probable) when the 
OPs are calculated from the {/max,i} sequence (from the 
{ATi} sequence). We also note that the patterns calcu¬ 
lated from the intensity peaks capture more determinism 
than those computed from the time intervals [notice the 
difference in the vertical scales of Figs. 2(a) and 2(b)]. 
This indicates that the timing of the high intensity peaks 
is more random than their peak values. 

The PE entropy. Fig. 2(c), quantifies this effect by 
decreasing sharply at the transition power. A similar 
behavior is observed when computing the HVG-entropy, 
Fig. 2(d). We have verified that, during the transition, 
the OP probabilities (both, for the intensity peaks and 
for the time intervals) are not consistent with the null hy¬ 
pothesis of full stochasticity: the probabilities lie outside 
the region consistent with random OPs, p ± 3ap, where 
p = 1/D! and dp = ^/p ( 1 — p )/N with N being the num¬ 
ber of data points In Fig. 2(c) the PE-entropy was 
computed from 0 = 3 OPs; a similar plot was obtained 
with 0 = 4 and 0 = 5 OPs (not shown). Larger O 


values were not considered due to the finite length of the 
dataset: while the “raw” intensity time series contains 50 
million data points, the number of high intensity peaks 
(above 2(t) shown in Fig. 1(d) is about 10^, depending 
on the pump power. 




FIG. 3. (a) The permutation entropy (left vertical axis) and 
the HVG-entropy (right vertical axis) vs the pump power. 
Galculations are performed for the ordinal patterns of length 
0 = 3. (b) Pdf-entropy calculated from the distribution of 
intensity values. 

The probabilities of patterns 012 and 210 provide a 
measure of the persistence of the time-series, i.e., the 
probability that the sign of Xi — Xi-i persists in the 
next step [l^. Thus, at the transition, if there are 
two consecutive peaks with increasing height, the next 
peak is likely to be larger than the previous one (and if 
there are two consecutive peaks with decreasing height, 
the next one is likely to be smaller than the previ¬ 
ous one); on the contrary, in the sequence of time- 
intervals, two consecutive intervals that are increasingly 
long (ATi < ATi+i) are likely to be followed by shorter 
interval (AT^+i > /S.Ti+ 2 ), and two consecutive decreas¬ 
ing intervals (AT^ > AT^+i) are likely to be followed by 
a longer one (AT^+i < ATi_|_ 2 ). 




FIG. 4. (a) Permutation entropy vs the lag-time before (0.85 
W), at (0.9 W), and after (0.95 W) the transition to optical 
turbulence, (b) Autocorrelation function vs lag, for the same 
pump powers as panel (a). 

Let us note that the optimal value of pattern length 
D depends of the length of correlations embedded into 
the time series @ . To investigate the presence of specific 
time-scales in the dynamics, we analyze the lagged time 
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series, i.e. the sequence of {/,,/i+T-,/i_|_ 2 r, ■ ■ ■}, where r 
is the lag time. We begin by considering the case r = 1 
(i.e., we analyze all data points). Figure 3(a) displays the 
PE entropy and the HVG entropy vs the pump power, 
and same behavior is seen in both entropies: there is a 
clear transition at pump power 0.9 W, where both en¬ 
tropies smoothly decrease. It is also observed that for 
the highest pump power, both entropies increase again. 
This reveals that during the transition there is an in¬ 
crease in the “ordering” of consecutive intensity values 
(that is captured by both entropies, which decrease), but 
for the highest pump power the trend reverses and the 
disorder increases. In contrast, the entropy computed in 
the conventional way and referred as pdf-entropy (i.e. the 
entropy calculated from the intensity pdfs of the initial 
intensity dynamics, I(t)) does not capture this behav¬ 
ior: as it can be seen in Fig. 3 (b), after the transition 
the pdf-entropy monotonously decreases with the pump 
power. Thus, SpE and Shvg provide consistent infor¬ 
mation, which complements that gained from the stan¬ 
dard pdf-entropy. The good agreement between the PE 
and HVG entropies, also seen in Eigs. 2(c) and 2(d), is 
remarkable because the two methods transform a time 
series into a sequence of integer numbers by using very 
different encoding rules. 

By varying the value of the lag time r, i.e. by taking 
into account not all points in data sets, but every second 
(r = 2) point, every third (r = 3) point etc., we are able 
to identify a specific oscillation time-scale in the intensity 
time-series during the transition. The PE entropy vs r 
for pump powers below (0.85 W), at (0.90 W) and above 
(0.95 W) the transition is displayed in Fig. 4(a). Here 
we can notice that, at the transition, there are specific 
lags for which the PE entropy decreases sharply. Similar 
results were obtained with the HVG entropy. 

The sharp minima indicate that, for pump power 0.90 
W and for specific lags, 6 different patterns of length 
D = 3 are not all equally probable, and thus, there 
are serial correlations in the sequence of lagged intensity 
values. To explore the length of such correlations, we 
computed the PE entropy using longer ordinary patterns 
{D = 4 and D = 6) and found that the minima were 
more pronounced, revealing the existence of long serial 
correlations. These correlations are not captured by the 
classical autocorrelation function shown in Fig. 4(b) for 
comparison purposes, that displays only a smooth varia¬ 
tion with T. 

To investigate the nature of these correlations we plot 
on Fig. 5(a) how the probabilities of 6 different ordinary 
patterns of length D = 3 depend on the lag time. We 
note a periodic alternation in which 012 and 210 became 
the more probable or the less probable patterns. The 
probabilities of the other four patterns are similar (no 
clear clusters are seen). The lag values for which 012 
and 210 are less probable correspond, as expected, to the 
lag values where the autocorrelation function is minimum 
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FIG. 5. (a) Probabilities of the ordinary patterns of length 
D = 3 vs lag time r. (b)-(d) Spatio-temporal structures iden¬ 
tified with the specific lags indicated with arrows in panel (a). 
The color scale indicates the value of I; with i = nr + j and 
r=396, 431 and 496 in units of the sampling time. The pump 
power is 0.9W. 


(and negative). However, unexpectedly, the lag values for 
which 012 and 210 are more probable, do not correspond 
to the maxima of the autocorrelation; and moreover, for 
the lag values where the autocorrelation is maxima, the 
all six ordinary patterns have similar probabilities. These 
observations suggest that ordinal analysis identifies sub¬ 
tle correlations in the ordering of data points, which are 
not seen by the standard autocorrelation function, that 
measures correlations in the values of the data points. 

These “order correlations” result into different types 
of spatio-temporal patterns. Let us recall that the ini¬ 
tial intensity spatio-temporal dynamics could be seen as 
intensity spatio-temporal dynamics 0- Similar concept 
has been recently employed for studying interaction of 
cavity solitons [1^ and topological solitons as address¬ 
able bits 0. Here we apply this concept and choose 
specific lag times defined from Fig.5(a) (shown by ar¬ 
rows) to be used as round-trip times in processing of the 
initial data series (see 0 for details of plotting spatio- 
temporal dynamics). Figures 5(b)-(d) display examples 
with lags such that patterns 012 and 210 are as probable 
[Fig. 5(b)], more probable [Fig. 5(c)]) and less probable 
[Fig. 5(d)] than the other four patterns. We can see that 
these spatio-temporal dynamics of patterns display clear 
and different coherent structures. These observations can 
be useful for confronting the predictions of state-of-the- 
art laser models with empirical data, and the theoretical 
studies could provide insight into the physical mecha- 
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nisms underlying these correlations. 

To summarize, by applying two independent tools of 
nonlinear time-series analysis we have uncovered long- 
range temporal correlations in the intensity output of 
a fibre laser during the transition to a wave turbulence 
regime. Output of laser radiation is easily measurable, 
making these easily implementable methods useful and 
valuable techniques for investigating coherent structures 
in complex laser radiation. Both approaches can be ap¬ 
plied to any high-dimensional complex systems that un¬ 
dergo similar transitions accompanied by the generation 
of extreme fluctuations. 
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